home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_5 / a5_4.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.1 KB  |  121 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 5.4 (Cubic Splines).
  9. % Section    5.3,    Interpolation by Spline Functions, Page 297
  10. echo on; clc; format short; hold off; clear
  11.  
  12. % This program finds the cubic spline interpolant.
  13.  
  14. % Given a set of data points
  15.  
  16. %      { (x , y ), (x , y ) ,..., (x , y ) }.
  17. %          1   1     2   2          n   n
  18.  
  19. % The abscissas and ordinates are stored in X and Y, respectively.
  20.  
  21. % X = [x , x  ,..., x ];  Y = [y , y  ,..., y ];
  22. %       1   2        n          1   2        n
  23.  
  24. % Remark. csfit.m  and cs.m  are used for Algorithm 5.4
  25.  
  26. pause % Press any key to continue.
  27.  
  28. clc;
  29. %    Place the abscissas for the points in  X.
  30.  
  31. %    Place the ordinates for the points in  Y.
  32.  
  33. X = [0   1     3     4     6  ];
  34.  
  35. Y = [0   0.5   2.0   1.5   2.5];
  36.  
  37. X    =    [0   1    2    3 ];
  38. Y    =    [0  0.5  2.0  1.5];
  39.  
  40. pause    % Press any key to graph data points.
  41.  
  42. clc; clg;
  43. a = min(X)-0.2; b = max(X)+0.2;
  44. c = min(Y)-0.2; d = max(Y)+0.2;
  45. axis([a b c d]);...
  46. plot(X,Y,'or');...
  47. hold on;...
  48. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  49. xlabel('x');...
  50. ylabel('y');...
  51. title('The given x-y data points');...
  52. grid;...
  53. axis;...
  54. hold off;...
  55. shg; pause    % Press any key to continue.
  56.  
  57. points = [X;Y];
  58. format short;
  59. clc,disp('The given x-y points:'),...
  60. disp('      x         y'),disp(points')
  61.  
  62. pause    % Press any key for the menu of spline curves.
  63.  
  64. clc; 
  65. Mx1='     Available curves are:';
  66. Mx2='(1)  Clamped Spline';
  67. Mx3='(2)  Natural Spline';
  68. Mx4='(3)  Extrapolated Spline (Cubic runout spline)';
  69. Mx5='(4)  Parabolically Terminated Spline';
  70. Mx6='(5)  Endpoint Curvature Adjusted Spline';
  71. clc,disp(''),disp(Mx1),disp(''),disp(Mx2),disp(''),disp(Mx3),,...
  72. disp(''),disp(Mx4),disp(''),disp(Mx5),disp(''),disp(Mx6),...
  73. ct = input('Enter the spline type <1-5> ');...
  74. if length(ct)==0, ct=2; end;...
  75. disp('');...
  76. S = csfit(X,Y,ct);
  77.  
  78. pause    % Press any key to fit the spline curve.
  79.  
  80. clc;
  81. hs = (max(X)-min(X))/200;
  82. Xs = min(X):hs:max(X);
  83. Ys = cs(S,X,Xs);
  84. axis([a b c d]);...
  85. plot(X,Y,'or',Xs,Ys,'-g');...
  86. hold on;...
  87. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  88. xlabel('x');...
  89. ylabel('y');...
  90. title('The cubic spline: y = S(x)');...
  91. grid;...
  92. axis;...
  93. hold off;...
  94. shg; pause    % Press any key to continue.
  95.  
  96. Mx1 ='The given x-y points:';
  97. Mx2 = '      x         y';
  98. if ct==1,Mx='clamped spline';end
  99. if ct==2,Mx='natural spline';end
  100. if ct==3,Mx='extrapolated spline (cubic runout spline)';end
  101. if ct==4,Mx='parabolically terminated spline';end
  102. if ct==5,Mx='endpoint curvature adjusted spline';end
  103. Mx3 = ['The ',Mx,' was determined.'];
  104. Mx4 = 'The spline coefficients are:'; 
  105. clc,echo off,diary output,...
  106. disp(Mx1),disp(Mx2),disp([X;Y]'),...
  107. disp(Mx3),disp(Mx4),disp(''),...
  108. [n1 n2] = size(S);...
  109. for k = 1:n1,
  110.   Mx5 = ['S',num2str(k),'(x):'];
  111.   disp(Mx5);disp(S(k,4:-1:1));
  112. end,diary off,echo on
  113.  
  114. pause % Press any key to continue.
  115.  
  116. points = [X;Y;CS(S,X,X);Y-CS(S,X,X)]';
  117. Mx5 = '    x(k)      y(k)      S(x(k))   error';
  118. clc,echo off,diary output,...
  119. disp(''),disp(Mx5),disp(points),diary off,echo on
  120.  
  121.